This software has been approved for release by the U.S. Geological Survey (USGS). Although the software has been subjected to rigorous review, the USGS reserves the right to update the software as needed pursuant to further analysis and review. No warranty, expressed or implied, is made by the USGS or the U.S. Government as to the functionality of the software and related material nor shall the fact of release constitute any such warranty. Furthermore, the software is released on condition that neither the USGS nor the U.S. Government shall be held liable for any damages resulting from its authorized or unauthorized use.
This vignette is intended to demonstrate how to run a multi-species Bayesian hierarchical model, developed at the U.S. Geological Survey, which predicts bat fatalities at U.S. wind facilities. We show how to load included data objects, run the model, and generate bat fatality predictions.
This repository is installed as an R package and has been primarily developed to support efficient data access and scientific reproducibility.
Several prepared data frames are included in the package in order to easily run this analysis.
The first data frame describes turbines recorded in the U.S. Wind Turbine
Database (USGS), including both those currently in operation and
those that have been decommissioned. Load the data frame and its
documentation using the following commands. The code to prepare these
data can be accessed at: data-raw/wind/turbines.R file.
#Load the turbines data frame
#The preview in this vignette shows the first 10 rows only
NationalBatsAndWind::turbines
#Run the following command to view the metadata for each column in `turbines`
#This will appear within the "Help" tab to the right in RStudio
?NationalBatsAndWind::turbines
The second data frame aggregates turbines by wind energy facilities.
The data frame includes descriptive statistics for turbine
specifications and other covariates associated with each facility.
Facilities are grouped by development phase and decommission status. For
example, if a facility has multiple phases, the facility may appear in
multiple rows, but the start_year and end_year
values for these phases will vary. This occurs because a single facility
may have constructed and/or decommissioned wind turbines across several
years. Load this data frame and its documentation using the following
commands. The code to prepare these data are in the
data-raw/wind/facilities.R file.
#Load the facilities data frame
#The preview in this vignette shows the first 10 rows only
NationalBatsAndWind::facilities
#Run the following command to view the metadata for each column in `facilities`
?NationalBatsAndWind::facilities
We processed spatial covariates on a 11.5 x 11.5 km grid across the
contiguous U.S. The original grid was created by the National Renewable
Energy Laboratory (NREL) to estimate and display wind energy development
potential under a variety of conservation, social, and geographic
restrictions. Processed covariate data across the grid is stored in
data/grids/nrel/final_nrel_grid_with_covariates.shp. Use
the following commands to load the processed grid and its
documentation.
#Load the nrel_grid data frame
#The preview in this vignette shows the first 10 rows only
NationalBatsAndWind::nrel_grid
#Run the following command to view the metadata for each column in `nrel_grid`
?NationalBatsAndWind::nrel_grid
To reduce the package size, the raw data layers to process these
covariates are not directly included in this repository. The code to
reproduce this spatial grid can be found at
data-raw/grid/add_covariates_to_nrel_grid.R. To rerun this
script fully, please reach out to the point of contact for this
repository to request copies of the raw data layers.
As an example of the data stored in the nrel_grid, we
will plot average windspeed at 10 meters above ground across the
continental U.S below.
#Load state data for outlines in map
state_lookup <- data.frame(state = state.abb, state_name = state.name)
states <- map_data("state")
states <- st_as_sf(states, coords = c('long', 'lat'), crs = 4326)
states <- st_transform(states, crs = st_crs(NationalBatsAndWind::nrel_grid))
states$x <- st_coordinates(states)[,1]
states$y <- st_coordinates(states)[,2]
states$state <- states$region
#Plot a continuous covariate from the NREL grid data
covariate_plot <- ggplot() +
#Select windspeed at 10 m above ground
geom_sf(
data = NationalBatsAndWind::nrel_grid,
aes(fill = windspeed_10_m_av),
color = "transparent") +
#Add state outlines
geom_polygon(
data = states,
aes(x = x, y = y, group = group),
color = "white", fill = "transparent", linewidth = 0.25) +
#Set graphical parameters for a pretty plot
scale_fill_gradientn(
colors = gist_ncar(100),
labels = scales::unit_format(unit = "m/s", big.mark = ","),
na.value="gray50") +
ggthemes::theme_map() +
labs(fill = "Average windspeed") +
theme_minimal() +
theme(plot.background = element_rect(fill = "white", color = "white"),
legend.position = "right",
text = element_text(size=16, colour="black"),
axis.text=element_text(size=12, colour="black"),
axis.title=element_text(size=12,face="bold")) +
ylab("") +
xlab("") +
coord_sf(crs = crs(NationalBatsAndWind::nrel_grid))
#Save plot to folder
ggsave(
"plots/nrel_grid_covariates_example.bmp",
plot = covariate_plot,
width = 9,
height = 4,
units = c("in"),
dpi = 300,
create.dir = T
)
#Print plot to vignette
covariate_plot
The bat fatality data in this software release originated from post-construction monitoring of bat fatalities at wind facilities across the United States. Because these data have been collected and delivered to the U.S. Fish & Wildlife Service (USFWS) related to the issuance of incidental take permits, there is often a requirement the data and results are available to the public. Each of the reports included in this release has previously been made publicly available, and we provide the public hyperlink for each report. The data in this release are formatted to support a reproducible analysis. Many other post-construction monitoring reports exist but are not included in this release due to business confidential rights/restrictions which may make them exempt from disclosure under Freedom of Information Act (FOIA) Exemption 4.
In our full analysis, which is not included here, we are able to run
our model across all reports shared with the USFWS and USGS, but in this
repository, we only include publicly available data. We plan to publish
a manuscript on an analysis of our full data set soon. For now, we are
releasing the software and a demo analysis on the subset of public data.
These public data are stored in the public_pcm data frame
and can be used to run a demo version of our model below. By showing how
these data are formatted, we hope others with additional PCM data that
we may not have access to will be able to prepare their PCM data in the
format shown in public_pcm and run our model on their own
data to independently verify our results. The scripts in
data-raw can help prepare additional PCM data in the
public_pcm format and here we relate covariate values from
the spatial grid to the PCM reports.
#Load mortality data
#The preview in this vignette shows the first 10 rows only
NationalBatsAndWind::public_pcm
#Check out the documentation for `public_pcm`
?NationalBatsAndWind::public_pcm
The following describes the model implemented in
inst/jags/multispecies_nmixture_without_ecoregion_poisson.txt.
The second version
(inst/jags/multispecies_nmixture_with_ecoregion_poisson.txt)
includes a random effect for ecoregion in the model. The
Poisson-Binomial mixture model for the count of carcasses \(Y_{fs}\) observed at facility \(f\) for species \(s\), out of \(N_{fs}\) total estimated fatalities, is
given by:
\[ \begin{aligned} Y_{fs} \ | \ N_{fs}, p_{f} &\overset{indep.}{\sim} \text{Binomial}\left(N_{fs}, p_{f}\right), \\ N_{fs} \ | \ \lambda_{fs} &\overset{indep.}{\sim} \text{Poisson}\left(\lambda_{fs}\right). \end{aligned} \]
A prior for the detection probability \(p_{f}\) is informed by the number of fatalities observed over the total bat fatalities originally estimated for facility \(f\) in its associated post-construction monitoring (PCM) report. Most of these PCM reports estimated total bats killed per megawatt, so we adjusted the estimates by the total rated capacity of the facility when applicable. These reports also used a variety of fatality rate estimators, some of which are known to be biased. Instead of using these potentially biased estimates as a fully trusted value, we borrowed the information to weakly inform the prior of the detection probability \(p_f\) at each facility. Parameters \(a_{f}\) and \(b_{f}\) were calculated to center the empirical detection probability of number of observed fatalities out of the total expected fatalities between the upper and lower uncertainty bounds calculated by the original estimator such that:
\[ p_{f} \sim \text{Beta}(a_{f}, b_{f}) \]
Then, based on the number of fatalities observed of each species at each facility and the weakly informed prior on the detection probability, the rate of fatalities \(\lambda_{fs}\) is modeled by a Bayesian GLMM with an offset based on the number of turbines at the facility \(n_{turbines}\):
\[ \begin{aligned} \log(\lambda_{fs}) = \text{log}(n_{turbines}) + \beta_s +& \\ & X_{fs,\text{seasons}}\beta_{\text{seasons}} + \\ & (1-I(\text{yes (1)/no (0) estimator adjusts for distance})) \times \frac{1}{\text{distance}_{f}} \beta_{\text{distance}} + \\ & X_{f, \text{covariates}}\beta_{\text{covariates,s}} + \\ & \epsilon_{\text{ecoregion}} + \epsilon_{\text{obs}} \end{aligned} \]
Each species has a species-specific fixed-effect intercept. Recall that JAGS parameterizes Normal distributions using mean and precision (a precision of 0.01 results in a very flat Normal distribution).
\[ \begin{aligned} \beta_{s} \sim \text{Normal}(0, 0.01) \end{aligned} \]
Each report has an adjustment for how much of the spring, summer, and fall seasons were searched. We assumed that if more of each season was searched, then more bat fatalities would be observed, so each \(\beta_{\text{seasons}}\) coefficient has a half-Normal prior with shared hyper-means and precision across species.
\[ \begin{aligned} \beta_{\text{season}_i, s} \sim \text{Normal}(\mu_{\text{season}_i}, \tau_{\text{season}}) T(0,) \end{aligned} \]
Some reports used an estimator that did not adjust for the number of bat carcasses that may have fallen outside of the search radius. For those reports, an adjustment is included in the model to account for the density of carcasses that may have fallen outside the searched area.
\[ \begin{aligned} (1-I(\text{yes (1)/no (0) estimator adjusts for distance})) \times \frac{1}{\text{distance}_{f}} \beta_{\text{distance}} \\ \end{aligned} \]
Species-specific coefficients for the covariate relationships \(X_{f, \text{covariates}}\beta_{\text{covariates,s}}\) were assumed to have covariate-specific hyper-means and precision.
\[ \begin{aligned} \beta_{\text{covariate}_j, s} \sim \text{Normal}(\mu_{\text{covariate}_j}, \tau_{\text{covariate}_j}) \end{aligned} \]
Additionally, there are optional random effects for ecoregion \(\epsilon_{\text{ecoregion}}\) and observation-level random effects \(\epsilon_{\text{obs}}\) to help account for overdispersion.
Other parameters in the model were assumed to have the following priors:
\[ \begin{aligned} \beta_{\text{distance}} &\sim \text{Normal}(0, 0.01) \\ \mu_{\text{season}_i} &\sim \text{Normal}(0, 0.01) T(0,) \\ \tau_{\text{season}} &\sim 1 / \sqrt{\text{Uniform}(0, 10)} \\ \mu_{\text{covariate}_j} &\sim \text{Normal}(0, 0.01) \\ \tau_{\text{covariate}_j} &\sim 1 / \sqrt{\text{Uniform}(0, 10)} \end{aligned} \]
There are a few additional steps for processing the data that are required before running the model. The steps that we will process in the next steps include:
In this step, we join the species ranges with each PCM study to determine whether the facility falls within each species range. This step will help the model estimate species-specific covariate effects only when a facility falls within the expected range of the species or fatalities of the species were observed even though the facility is out of the known range.
#Load shapefile of species ranges
spc_ranges <- read_sf("../data-raw/bats/ranges/USGS_ScienceBase_coded.shp") %>%
st_transform(crs = st_crs(NationalBatsAndWind::nrel_grid))
#Determine species codes that exist in the public_pcm data
spc_start <- which(colnames(public_pcm) == "EPFU (Big brown)")
spc_end <- which(colnames(public_pcm) == "NYFE (pocketed fre-tailed)")
spcs <- colnames(public_pcm)[spc_start:spc_end]
spcs_code <- sapply(spcs, substr, 1, 4)
#Select the species ranges that exist within the public_pcm data
spc_ranges <- spc_ranges %>%
filter(Spp_code %in% spcs_code) %>%
st_transform(crs = st_crs(NationalBatsAndWind::nrel_grid))
#Add species ranges to the mortality data (1 if facility in range, 0 otherwise)
public_pcm_with_spc <- NationalBatsAndWind::public_pcm |>
sf::st_as_sf(coords = c("x", "y"), crs = st_crs(NationalBatsAndWind::nrel_grid)) |>
sf::st_join(spc_ranges) |>
dplyr::select(-SCI_NAME, -COMMON_NAM) |>
tidyr::pivot_wider(
names_from = Spp_code,
values_from = Spp_code,
values_fn = ~1,
values_fill = 0) |>
st_drop_geometry() |>
#Join nrel_grid with the facility data
dplyr::left_join(NationalBatsAndWind::nrel_grid, by = "fid")
Now we prepare the covariate data matrix for the facilities with PCM
reports. In summary, a few of the covariates in the
public_pcm are slightly altered or renamed to be consistent
with the prediction steps later on. Additionally, several other
covariate vectors that help with the bias adjustments for search
strategies are created for the model.
#Create the covariate matrix
Xdf <- public_pcm_with_spc |>
st_drop_geometry() |>
#Edit values of a few covariates to account for zeroes and missing values
dplyr::mutate(
years_since_wns = uswtdb_min_p_year - wns_yoa_av,
years_since_wns_nonneg = if_else(years_since_wns < 0, 0, years_since_wns),
max_cutin_speed = suppressWarnings(as.numeric(max_cutin_speed)),
adj_distance = grepl("no", estimator_accounts_for_max_dist) * 1,
max_cutin_speed = if_else(is.na(max_cutin_speed), 3, max_cutin_speed),
dist_searched_transform = 1/max_distance_searched) |>
#Select covariates to include in model matrix
dplyr::select(
#Select covariates that are not scaled first
fid, ecoregion2_mode, used_cutin,
adj_distance, dist_searched_transform,
prop_spring, prop_summer, prop_fall,
n_turbine = estimate_relevant_tnum,
#Rename and select these covariates
total_rsa = estimate_relevant_rsa,
total_mw = estimate_relevant_cap,
#Also select the following
starts_with("nlcd_ilr"),
mean_rsa, mean_lrs, age,
lrs_mean_windspeed_av, physdiv_av,
vapor_pr_av
)
#Some estimators account for fatalities that may fall outside of the search area
#This is a boolean vector that indicates which reports still need adjusting
adj_distance <- Xdf$adj_distance
Here, we center and scale most of the covariates and save the parameter values for centering and scaling to use again during prediction.
#Center and scale the covariate matrix
Xdfcs <- Xdf |>
#Center and scale all variables except those that are factors or indicators
mutate(across(-c(1:9), \(x) scale(x, center = TRUE, scale = sd(x, na.rm = TRUE)))) |>
#The season proportions don't handle scaling well. So center to 1 so that
#a full season search drops out of the model (this is what we want the
#model to estimate by default) and "scale" by 1 so no-scaling occurs
mutate(across(starts_with("prop"), \(x) scale(x, center = 1, scale = 1))) |>
#The distance_searched_transform has a half-normal prior, so let's set
#the center to 0 (since estimates that don't need to be distance adjusted
#will have a value of 0 and drop out of the model anyway) and scale
mutate(across(dist_searched_transform, \(x) scale(x, center = 0, scale = sd(x, na.rm = TRUE)))) |>
#Calculate quadratic effects for some variables after center/scale operation
mutate(vapor_pr_av_sq = vapor_pr_av^2)
#Save the center and scale information for the covariate matrix to use later
cs_save <- data.frame(rbind(
c('used_cutin', 0, 1),
unlist(c('dist_searched_transform', attributes(Xdfcs$dist_searched_transform)[2:3])),
unlist(c('prop_spring', attributes(Xdfcs$prop_spring)[2:3])),
unlist(c('prop_summer', attributes(Xdfcs$prop_summer)[2:3])),
unlist(c('prop_fall', attributes(Xdfcs$prop_fall)[2:3])),
unlist(c('total_rsa',attributes(Xdfcs$total_rsa)[2:3])),
unlist(c('mean_rsa',attributes(Xdfcs$mean_rsa)[2:3])),
unlist(c('mean_lrs',attributes(Xdfcs$mean_lrs)[2:3])),
unlist(c('age',attributes(Xdfcs$age)[2:3])),
unlist(c('lrs_mean_windspeed_av', attributes(Xdfcs$lrs_mean_windspeed_av)[2:3])),
unlist(c('physdiv_av',attributes(Xdfcs$physdiv_av)[2:3])),
unlist(c('vapor_pr_av',attributes(Xdfcs$vapor_pr_av)[2:3])),
unlist(c('vapor_pr_av_sq',attributes(Xdfcs$vapor_pr_av)[2:3])),
unlist(c('nlcd_ilr_1',attributes(Xdfcs$nlcd_ilr_1)[2:3])),
unlist(c('nlcd_ilr_2',attributes(Xdfcs$nlcd_ilr_2)[2:3])),
unlist(c('nlcd_ilr_3',attributes(Xdfcs$nlcd_ilr_3)[2:3])),
unlist(c('nlcd_ilr_4',attributes(Xdfcs$nlcd_ilr_4)[2:3])),
unlist(c('nlcd_ilr_5',attributes(Xdfcs$nlcd_ilr_5)[2:3])),
unlist(c('nlcd_ilr_6',attributes(Xdfcs$nlcd_ilr_6)[2:3])),
unlist(c('nlcd_ilr_7',attributes(Xdfcs$nlcd_ilr_7)[2:3])),
unlist(c('nlcd_ilr_8',attributes(Xdfcs$nlcd_ilr_8)[2:3])),
unlist(c('nlcd_ilr_9',attributes(Xdfcs$nlcd_ilr_9)[2:3])))) |>
mutate(
scaled.center = as.numeric(scaled.center),
scaled.scale = as.numeric(scaled.scale)) |>
rename(colname = V1)
#Update rownames of the saved data and remove the names in the first column
rownames(cs_save) <- cs_save$colname
cs_save <- cs_save[,2:3]
#Some reports require adjustments for total amount of each season searched
S <- Xdfcs |> dplyr::select(prop_spring, prop_summer, prop_fall)
#A vector of the max distance searched from turbines in each report
dist_searched_transform <- as.numeric(Xdfcs$dist_searched_transform)
Using the centered and scaled covariate values, we set up the model.matrix of covariates for the model by selecting those to include. We also select the species to include in the model and subset the species ranges based on the selected species.
#Use an offset for the number of turbines relevant to each estimate
offset <- Xdf$n_turbine
#Set up the model matrix by selecting covariates to include in the analysis
X <- model.matrix(
~ 1 +
mean_rsa +
age +
physdiv_av +
mean_lrs +
lrs_mean_windspeed_av +
nlcd_ilr_1 +
nlcd_ilr_2 +
nlcd_ilr_3 +
nlcd_ilr_4 +
nlcd_ilr_5 +
nlcd_ilr_6 +
nlcd_ilr_7 +
nlcd_ilr_8 +
nlcd_ilr_9,
Xdfcs)
#Remove the intercept from the model matrix, because a new intercept coefficient
#will be estimated for each species.
if(ncol(X) == 2){
X <- as.matrix(X[,2], nrow = nrow(public_pcm_with_spc), ncol = 1)
} else {
X <- X[,2:ncol(X)]
}
#Select the species to include in the model
species_columns <- c(
"LACI (Hoary)",
"LABO (Eastern red)",
"EPFU (Big brown)"
)
#Subset the species range columns based on the selected species
R <- public_pcm_with_spc[,substr(species_columns, 1, 4)] |>
st_drop_geometry()
The data objects that were created above can be used as inputs to the multi-species model that we described in Model Details. The function-call below is set to a low number of iterations in order to demonstrate how the model runs in a quick fashion. However, in order to achieve convergence, a much higher number of iterations is suggested. In the comments, we provide a suggested number of iterations that we have used with this model to reach approximate convergence. Running the code below should take approximately 10 minutes on a standard 2020s computer. Expect the model to require up to several hours when increasing the number of iterations to achieve convergence.
Please keep in mind when running this model and reviewing the analysis, that the model is run only on the publicly available data included in this release. In the near future we plan to publish our results on the full data set, which includes business confidential reports. Our analyses have shown that the results on the full data set are more robust and don’t necessarily agree fully with the example results described below.
set.seed(20250404)
#Begin the model run with the processed inputs
#We recommend running the model for more iterations than specified in the demo
#Consider replacing the jags_pars with the following parameters
#These parameters will produce better convergence but require
#more runtime (several hours)
# jags_pars = list(
# ni = 1000000,
# nc = 3,
# na = 50000,
# nb = 50000,
# nt = 10,
# parallel = TRUE),
system.time({
jl <- NationalBatsAndWind::run_multispecies_model(
public_pcm_with_spc,
species_columns,
X,
R,
S,
offset,
adj_distance,
dist_searched_transform,
cs_save,
jags_pars = list(
ni = 100000,
nc = 3,
na = 50000,
nb = 50000,
nt = 3,
parallel=TRUE),
ecor = FALSE,
model = "nmixture",
family = "poisson",
params = c(
"beta_range", "beta_season", "beta_dist", "beta", "beta_mu",
"eps", "eps_sigma"
))})
#Save the JAGS outputs in two objects
jags_data <- jl$jags_data
jags_fit <- jl$jags_fit
These outputs can also be saved to a file to load at another time.
#If desired, save the outputs to a file
save(jags_data, file = "../data/models/model_run_data.rda")
save(jags_fit, file = "../data/models/model_run_fit.rda")
Previous model outputs can be reloaded. If you are returning to this vignette, please run all code chunks prior to fitting and saving the model (the previous two code chunks), so that all R objects prior to the model run are available.
#If desired, load previous model outputs without having to rerun the models
#Make sure to run the other code chunks before the "Run the model" section also
load("../data/models/model_run_data.rda")
load("../data/models/model_run_fit.rda")
Traceplots and whiskerplots can help evaluate convergence on model parameters. For a model run with a low number of iterations and a high number of parameters, it isn’t likely that convergence will be reached. Traceplots for each parameter should look like a “field of grass” that is centered around one mean across all iterations. Another sign of convergence is R-hat values < 1.1 or ideally < 1.01.
#Consider checking the traceplots of the following parameters for convergence.
#Edit the parameters variable in the traceplot() function below:
#parameters="beta_range"
#parameters="beta_dist"
#parameters="beta_season"
#parameters="beta"
jagsUI::traceplot(jags_fit, parameters = "beta_range")
Our multi-species model of bat fatalities at wind energy facilities estimates a variety of parameters that can be summed together to predict species-specific fatality rates. First, we estimate the average annual rate of bats killed per turbine for each species. In the model, this parameter is \(\beta_s\). This parameter can be interpreted as the number of bats killed per turbine after marginalizing across the mean of all other covariate effects.
#Use the simulated parameter values to determine the mean and lower/upper bounds
#of the 95% credible interval for each intercept parameter
range_intercepts_data <- data.frame(
covariates = species_columns,
color = jags_fit$overlap0$beta_range,
lo = apply(exp(jags_fit$sims.list$beta_range), 2, quantile, 0.025),
mean = apply(exp(jags_fit$sims.list$beta_range), 2, mean),
up = apply(exp(jags_fit$sims.list$beta_range), 2, quantile, 0.975))
#Create a pretty plot of the intercept coefficients
range_intercepts_plot <- ggplot(range_intercepts_data) +
geom_errorbar(aes(x = covariates, ymin = lo, ymax = up)) +
geom_point(aes(x = covariates, y = mean)) +
geom_hline(yintercept = 0) +
ylim(min(c(range_intercepts_data$lo, 0)), max(range_intercepts_data$up)) +
coord_flip() +
theme_minimal() +
#ggtitle("Average annual fatalities per turbine") +
theme(plot.background = element_rect(fill = "white", color = "white"),
legend.position = "bottomleft",
text = element_text(size=16, colour="black"),
axis.text=element_text(size=12, colour="black"),
axis.title=element_text(size=12,face="bold")) +
ylab("Average annual fatalities per turbine") +
xlab("Species")
#Save the plot for later
ggsave(
"plots/multispecies_all_range_intercepts.png",
plot = range_intercepts_plot,
width = 6,
height = 6,
units = c("in"),
dpi = 300,
create.dir = T
)
#Print out the plot below
range_intercepts_plot
We can also estimate the average effect of each covariate in the model across all species. These parameters are considered hyper-means and correspond to the \(\mu_{\text{covariate}_j}\) parameters in the model.
#Use the simulated parameter values to determine the mean and lower/upper bounds
#of the 95% credible interval for each covariate effect hypermean
covariate_data <- data.frame(
covariates = colnames(jags_data$X),
color = jags_fit$overlap0$beta_mu,
lo = apply(jags_fit$sims.list$beta_mu[,], 2, quantile, 0.025),
mean = apply(jags_fit$sims.list$beta_mu[,], 2, mean),
up = apply(jags_fit$sims.list$beta_mu[,], 2, quantile, 0.975)) |>
mutate(Covariate = case_when(
covariates == "mean_rsa" ~ "Mean rotor swept area",
covariates == "age" ~ "Years since built",
covariates == "mean_lrs" ~ "Mean lower rotor sweep",
covariates == "physdiv_av" ~ "Physiographic diversity",
covariates == "lrs_mean_windspeed_av" ~ "Mean windspeed",
covariates == "nlcd_ilr_1" ~ "NLCD Ratio 1",
covariates == "nlcd_ilr_2" ~ "NLCD Ratio 2",
covariates == "nlcd_ilr_3" ~ "NLCD Ratio 3",
covariates == "nlcd_ilr_4" ~ "NLCD Ratio 4",
covariates == "nlcd_ilr_5" ~ "NLCD Ratio 5",
covariates == "nlcd_ilr_6" ~ "NLCD Ratio 6",
covariates == "nlcd_ilr_7" ~ "NLCD Ratio 7",
covariates == "nlcd_ilr_8" ~ "NLCD Ratio 8",
covariates == "nlcd_ilr_9" ~ "NLCD Ratio 9"
))
#Set the order of the covariates to appear in the plot
order <- c("Mean rotor swept area",
"Years since built",
"Mean lower rotor sweep",
"Physiographic diversity",
"Mean windspeed",
"NLCD Ratio 1",
"NLCD Ratio 2",
"NLCD Ratio 3",
"NLCD Ratio 4",
"NLCD Ratio 5",
"NLCD Ratio 6",
"NLCD Ratio 7",
"NLCD Ratio 8",
"NLCD Ratio 9")
#Create a pretty plot of the covariate effect coefficients
covariate_plot <- ggplot(covariate_data |> filter(Covariate %in% order)) +
geom_errorbar(aes(x = Covariate, ymin = lo, ymax = up, color = color)) +
scale_color_manual(values=c("FALSE"="red","TRUE"="blue")) +
geom_point(aes(x = Covariate, y = mean)) +
geom_hline(yintercept = 0) +
scale_x_discrete(limits = rev(order)) +
coord_flip() +
ylim(min(covariate_data$lo), max(covariate_data$up)) +
labs(color = "95% credible interval\noverlaps 0") +
theme_minimal() +
ggtitle("Effects over all species") +
theme(title = element_text(size = 12, colour = "black", ),
plot.background = element_rect(fill = "white", color = "white"),
legend.position = "bottomleft",
text = element_text(size=16, colour="black"),
axis.text=element_text(size=12, colour="black"),
axis.title=element_text(size=12,face="bold")) +
ylab("Standardized effect size") +
xlab("Covariate")
#Save the plot for later
ggsave(
paste0("plots/covariate_effects_all_species.png"),
plot = covariate_plot,
width = 6,
height = 6,
units = c("in"),
dpi = 300,
create.dir = T
)
#Print out the plot below
covariate_plot
The hyper-means do not reveal much information. However, in the next section, we will consider how these parameters could have strong species-level effects.
Let’s recreate the plot of covariate effects above, but plot how the effects vary for each species. These parameters are the species-specific effects in the model \(\beta_{\text{covariate}_j}\).
#Determine the species name (four letter code) from the species_columns
species_names <- substr(species_columns, 1, 4)
#For each selected species, create the plot
for(species in seq(1,length(species_names))) {
#Select the species
species_name <- species_names[species]
#Determine the range of the species
species_range <- spc_ranges |> filter(Spp_code == species_name)
#Use the simulated parameter values to determine the mean and bounds
#of the 95% credible interval for each species-specific covariate effect
covariate_data <- data.frame(
covariates = colnames(jags_data$X),
color = jags_fit$overlap0$beta[,species],
lo = apply(jags_fit$sims.list$beta[,,species], 2, quantile, 0.025),
mean = apply(jags_fit$sims.list$beta[,,species], 2, mean),
up = apply(jags_fit$sims.list$beta[,,species], 2, quantile, 0.975)) |>
mutate(Covariate = case_when(
covariates == "mean_rsa" ~ "Mean rotor swept area",
covariates == "age" ~ "Years since built",
covariates == "mean_lrs" ~ "Mean lower rotor sweep",
covariates == "physdiv_av" ~ "Physiographic diversity",
covariates == "lrs_mean_windspeed_av" ~ "Mean windspeed",
covariates == "nlcd_ilr_1" ~ "NLCD Ratio 1",
covariates == "nlcd_ilr_2" ~ "NLCD Ratio 2",
covariates == "nlcd_ilr_3" ~ "NLCD Ratio 3",
covariates == "nlcd_ilr_4" ~ "NLCD Ratio 4",
covariates == "nlcd_ilr_5" ~ "NLCD Ratio 5",
covariates == "nlcd_ilr_6" ~ "NLCD Ratio 6",
covariates == "nlcd_ilr_7" ~ "NLCD Ratio 7",
covariates == "nlcd_ilr_8" ~ "NLCD Ratio 8",
covariates == "nlcd_ilr_9" ~ "NLCD Ratio 9"
))
#Set the order of the covariates to appear in the plot
order <- c("Mean rotor swept area",
"Years since built",
"Mean lower rotor sweep",
"Physiographic diversity",
"Mean windspeed",
"NLCD Ratio 1",
"NLCD Ratio 2",
"NLCD Ratio 3",
"NLCD Ratio 4",
"NLCD Ratio 5",
"NLCD Ratio 6",
"NLCD Ratio 7",
"NLCD Ratio 8",
"NLCD Ratio 9")
#Create a pretty plot of the covariate effect coefficients
spc_covariate_plot <- ggplot(covariate_data |> filter(Covariate %in% order)) +
geom_errorbar(aes(x = Covariate, ymin = lo, ymax = up, color = color)) +
scale_color_manual(values=c("FALSE"="red","TRUE"="blue")) +
geom_point(aes(x = Covariate, y = mean)) +
geom_hline(yintercept = 0) +
scale_x_discrete(limits = rev(order)) +
coord_flip() +
ylim(min(covariate_data$lo), max(covariate_data$up)) +
labs(color = "95% credible interval\noverlaps 0") +
theme_minimal() +
ggtitle(paste0(species_name, "-specific effects")) +
theme(title = element_text(size = 12, colour = "black", ),
plot.background = element_rect(fill = "white", color = "white"),
legend.position = "bottomleft",
text = element_text(size=16, colour="black"),
axis.text=element_text(size=12, colour="black"),
axis.title=element_text(size=12,face="bold")) +
ylab("Standardized effect size") +
xlab("Covariate")
#Save the plot for later
ggsave(
paste0("plots/", species_name, "_all_covariate_effects.png"),
plot = spc_covariate_plot,
width = 6,
height = 6,
units = c("in"),
dpi = 300,
create.dir = T
)
}
The plots that we generate in the code above are saved to the
vignettes/plots folder. For each species, the code
generates a “????_all_covariate_effects.png” file, where the question
marks are replaced with each four-letter species code. Let’s look at the
species-specific covariate effects for hoary bat (LACI = Lasiurus
cinereus):
The species-specific results for hoary bat reveal stronger covariate effects. As the age of the facility and turbine ground clearance (mean lower rotor sweep) increase, there is a lower rate of hoary bat fatalities observed. Keep in mind, these are correlative observations, not necessarily causal and may be confounded with other covariates not included in the model. Additionally, there are a few landcover effects that are strongly associated with increased or decreased hoary bat fatality rates. These landcover effects are difficult to interpret directly, however, because they are based on a isometric logratio transformation which incorporates multiple landcovers into each covariate. In addition, please recall that this model evaluates only the publicly available data included in this release. In the near future we plan to publish our results on the full data set, which includes business confidential reports that we cannot share as raw data. Our analyses have shown that the results on the full data set are more robust and don’t necessarily agree fully with the example results described below. Finally, there may be some stochasticity in model runs, especially with a low number of iterations, so the strength of the effects described here may differ slightly.
Check out the additional
????_all_covariate_effects.png
Using the estimated covariate relationships, we can estimate total
bat fatalities for each species across a retrospective prediction for
all wind turbines that have been on the landscape from 2001 to 2024.
Recall that we compiled all turbines and wind facilities from the U.S.
Wind Turbine Database into the
NationalBatsAndWind::turbines and
NationalBatsAndWind::facilities data frames in this
package. For each of these facilities, we calculate the covariates used
in the model above, then use the covariate effects to predict the total
number of fatalities for each species that occurred at specific
facilities over each year.
#Determine the species name (four letter code) from the species_columns
species_names <- substr(species_columns, 1, 4)
#For each selected species, create the plot
for(species in seq(1,length(species_names))) {
#Select the species
species_name <- species_names[species]
#Determine the range of the species
species_range <- spc_ranges |> filter(Spp_code == species_name)
#Retrieve the center/scale information for the covariates
cs_save <- jags_data$cs_save
#Create a data.frame of parameter values to run predictions on
parameter_grid <- data.frame(intercept = rep(1, 24),
prop_spring = rep(1, 24),
prop_summer = rep(1, 24),
prop_fall = rep(1, 24),
cutin = rep(0, 24),
max_cutin_speed = rep(3, 24),
year = c(2001:2024))
#Edit sf settings to avoid warnings/errors
sf::sf_use_s2(FALSE)
#Determine which facilities fall within the selected species range
facilities_st_in_range <- st_join(facilities, species_range) |>
filter(!is.na(Spp_code)) |>
sf::st_drop_geometry() |>
left_join(NationalBatsAndWind::nrel_grid, "fid")
#Make of copy of the nrel_grid to save outputs to
nrel_grid_out <- NationalBatsAndWind::nrel_grid
#Print some messages to provide progress feedback to user
print(paste0("Processing species: ", species_name))
#For each row in the parameter grid
for(jj in 1:nrow(parameter_grid)){
#Filter to facilities in the species range to those that existed that year
facilities_st_filtered <- facilities_st_in_range |>
filter(start_year <= parameter_grid$year[jj] & end_year > parameter_grid$year[jj])
#Create covariate data to predict fatalities at each of these facilities
facilities_for_pred <- facilities_st_filtered |>
mutate(
year = parameter_grid$year[jj],
age = parameter_grid$year[jj] - start_year,
n_turbine = uswtdb_p_tnum,
p_cap = interp_p_cap,
sum_t_cap = interp_p_cap,
total_rsa = uswtdb_p_tnum * 3.14 * (interp_rd / 2)^2,
total_mw = interp_p_cap,
years_since_wns = parameter_grid$year[jj] - wns_yoa_av,
years_since_wns_nonneg = if_else(years_since_wns <= 0, 0, years_since_wns),
intercept = parameter_grid[jj,1],
prop_spring = parameter_grid[jj,2],
prop_summer = parameter_grid[jj,3],
prop_fall = parameter_grid[jj,4],
mean_rsa = 3.14 * (interp_rd / 2)^2,
mean_hh = interp_hh,
mean_lrs = interp_lrs,
vapor_pr_av_sq = NA,
cutin = parameter_grid[jj,5],
max_cutin_speed = parameter_grid[jj,6])
#Center and scale the corresponding variables that are in the model
for(i in which(rownames(jags_data$cs_save) %in% colnames(jags_data$X))){
#Drop units on the column if present
if(class(facilities_for_pred[[rownames(jags_data$cs_save)[i]]])=="units"){
facilities_for_pred[[rownames(jags_data$cs_save)[i]]] <- units::drop_units(facilities_for_pred[[rownames(jags_data$cs_save)[i]]])
}
#Center and scale the columns according to matches in cs_save
facilities_for_pred <- facilities_for_pred %>%
mutate(!!rownames(jags_data$cs_save)[i] := (!!sym(rownames(jags_data$cs_save)[i]) - jags_data$cs_save[i,1] )/jags_data$cs_save[i,2])
}
#After center and scaling calculate the quadratic covariates
facilities_for_pred <- facilities_for_pred |>
mutate(vapor_pr_av_sq = vapor_pr_av^2)
#Predict fatalities for each facility in facilities_for_pred
temp <- facilities_for_pred %>%
NationalBatsAndWind::multispecies_predict(
jags_data,
jags_fit,
species = species,
resimulate = F,
fast = F)
predno_q2.5 <- paste0('pred_q2.5_', parameter_grid$year[jj])
predno_q50 <- paste0('pred_q50_', parameter_grid$year[jj])
predno_mean <- paste0('pred_mean_', parameter_grid$year[jj])
predno_q97.5 <- paste0('pred_q97.5_', parameter_grid$year[jj])
#Calculate the 0.025, 0.5, and 0.975 quantiles and mean of the posteriors
#across all fid cells included in the species range from the nrel_grid
out <- data.frame(fid = temp$cells_modeled)
out <- out %>%
mutate(q2.5 = apply(temp$pred, 1, quantile, 0.025)) %>%
mutate(q50 = apply(temp$pred, 1, quantile, 0.5)) %>%
mutate(mean = apply(temp$pred, 1, mean)) %>%
mutate(q97.5 = apply(temp$pred, 1, quantile, 0.975)) |>
group_by(fid) |>
#If there are multiple facilities in one fid, add these values together
summarize(!!predno_q2.5 := sum(q2.5),
!!predno_q50 := sum(q50),
!!predno_mean := sum(mean),
!!predno_q97.5 := sum(q97.5))
#Add the predictions to the copy of the nrel_grid
nrel_grid_out <- nrel_grid_out |>
left_join(out, by = 'fid')
}
#Calculate cumulative posterior medians for each state across all years
state_totals <- data.frame(nrel_grid_out) |>
group_by(state) |>
summarize(across(starts_with("pred_q50"), \(x) sum(x, na.rm = TRUE))) |>
mutate(state = tolower(state))
#Calculate lower bound of 95% credible interval for each year
year_totals_q2.5 <- data.frame(nrel_grid_out) |>
summarize(across(starts_with("pred_q2.5"), \(x) sum(x, na.rm = TRUE)))
#Calculate median for each year
year_totals_q50 <- data.frame(nrel_grid_out) |>
summarize(across(starts_with("pred_q50"), \(x) sum(x, na.rm = TRUE)))
#Calculate upper bound of 95% credible interval for each year
year_totals_q97.5 <- data.frame(nrel_grid_out) |>
summarize(across(starts_with("pred_q97.5"), \(x) sum(x, na.rm = TRUE)))
#Prepare a table of annual lower bounds, medians, and upper bounds
totals <- cbind(
c("lower", "median", "upper"),
rbind(as.vector(year_totals_q2.5[1,]),
as.vector(year_totals_q50[1,]),
as.vector(year_totals_q97.5[1,])))
colnames(totals) <- c("year", 2001:2024)
totals <- t(totals)
colnames(totals) <- totals[1,]
totals <- totals[-1,]
totals <- data.frame(totals)
totals$year <- 2001:2024
totals$lower <- as.numeric(totals$lower)
totals$median <- as.numeric(totals$median)
totals$upper <- as.numeric(totals$upper)
totals$year <- as.numeric(totals$year)
#Create an empty list of plots for storage
plots <- list()
#For each selected year create a plot of total fatalities per state
for(year in c(2001,2005,2010,2015,2020,2024)) {
filter_year <- year
#Select the facilities that exist on the landscape during filter_year
facilities_to_plot <- facilities |>
filter(start_year <= filter_year & end_year > filter_year) |>
st_transform(crs = 4326)
#Create a pretty plot of the total number of bat fatalities per state
#Include state outlines and points for each wind facility
output <- states |>
left_join(state_totals, by = "state") |>
ggplot() +
geom_polygon(aes(
x = x,
y = y,
group = group,
fill = !!sym(paste("pred_q50", year, sep = "_"))),
color = "white",
linewidth = 0.25) +
ggthemes::theme_map() +
labs(fill = "Total bats killed") +
ggtitle(paste("Year: ", year, sep = "")) +
scale_fill_gradientn(
colours=c("navyblue", "darkmagenta", "darkorange1", "gold"),
limits = c(1,max(state_totals |> dplyr::select(starts_with("pred_q50")))*1.1),
labels = scales::unit_format(unit = "", big.mark = ","),
na.value="gray50") +
ggthemes::theme_map() +
labs(fill = "Bat fatalities") +
theme_minimal() +
theme(
plot.background = element_rect(fill = "white", color = "white"),
legend.position = "right",
text = element_text(size=16, colour="black"),
axis.text=element_text(size=12, colour="black"),
axis.title=element_text(size=12,face="bold")) +
geom_point(
data = data.frame(cbind(x = c(-100), y = c(40), total = c(NA))),
aes(x = x, y = y, size = total, shape = "Wind facility"),
colour = "gray75") +
geom_sf(
data = facilities_to_plot,
aes(shape = "Location", color="white"),
color = "white", shape = 1, size = 0.25, show.legend = "point") +
guides(
shape = guide_legend(
"",
override.aes = list(shape = 1, size = 2, color = "black"))) +
ylab("") +
xlab("") +
coord_sf(crs = crs(NationalBatsAndWind::nrel_grid))
#Add plot to list of plots to use with gifski
plots[[length(plots)+1]] <- output
#Save plot for later
ggsave(
paste0("plots/", species_name, "_state_fatalities_", year, ".png"),
plot = output,
width = 9,
height = 4,
units = c("in"),
dpi = 300,
create.dir = T
)
}
#Combine plots into a .gif
library(gifski)
png_files <- list.files(
"plots/",
pattern = paste0(species_name, "_state_fatalities_.*png$"),
full.names = TRUE)
gifski(
png_files,
gif_file = paste0(
"plots/",
species_name,
"_state_fatalities_animation.gif"),
width = 2700,
height = 1200,
delay = 1)
#Calculate the cumulative sum of each statistic across years
totals_cumulative <- totals |>
mutate(
lower_cumul = cumsum(lower),
median_cumul = cumsum(median),
upper_cumul = cumsum(upper))
#Output the statistics to a file
write.csv(
totals_cumulative,
file = paste0("plots/", species_name, "_totals_cumulative.csv"))
#Create a ribbon plot of total annual fatalities
annual <- ggplot(totals, aes(x=year, y=median)) +
geom_point() +
geom_line() +
geom_ribbon(aes(ymin=lower, ymax=upper), linetype=2, alpha=0.1) +
scale_y_continuous(
limits = c(1,max(totals$upper)*1.1),
labels = scales::unit_format(unit = "", big.mark = ","),
na.value="gray50") +
theme_minimal() +
xlab("Year") +
ylab("Annual bat fatalities") +
theme(plot.background = element_rect(fill = "white", color = "white"),
legend.position = "right",
text = element_text(size=16, colour="black"),
axis.text=element_text(size=12, colour="black"),
axis.title=element_text(size=12,face="bold"))
#Save the plot for later
ggsave(
paste0("plots/", species_name, "_annual_fatalities.png"),
plot = annual,
width = 9,
height = 4,
units = c("in"),
dpi = 300,
create.dir = T
)
#Create a plot of total cumulative fatalities over time
cumulative <- ggplot(totals_cumulative, aes(x = year, y = median_cumul)) +
geom_point() +
geom_line() +
geom_ribbon(
aes(ymin = lower_cumul, ymax = upper_cumul),
linetype = 2, alpha = 0.1) +
scale_y_continuous(
limits = c(1, max(totals_cumulative$upper_cumul) * 1.1),
labels = scales::unit_format(unit = "", big.mark = ","),
na.value="gray50") +
theme_minimal() +
xlab("Year") +
ylab("Cumulative bat fatalities") +
theme(
plot.background = element_rect(fill = "white", color = "white"),
legend.position = "right",
text = element_text(size=16, colour="black"),
axis.text=element_text(size=12, colour="black"),
axis.title=element_text(size=12,face="bold"))
#Save plot for later
ggsave(
paste0("plots/", species_name, "_cumulative_fatalities.png"),
plot = cumulative,
width = 9,
height = 4,
units = c("in"),
dpi = 300,
create.dir = T
)
}
Let’s again consider just the hoary bat (LACI). But consider
reviewing the figures for the other species in the
vignettes/plots/ folder.
First, we can plot the estimated number of hoary bat fatalities per state over time. The animation shows how annual fatalities may change over time. Our predictions become more refined with our full data set and more covariates.
The model predicts that total fatalities have reached several hundred thousand hoary bats annually in recent years (95% credible interval).
Cumulatively, this example version of the model estimates that several million hoary bats have been killed due to wind energy nationwide since 2001.
Running the code chucks above will produce additional plots for other
bat species, including the big brown bat (EPFU = Eptesicus
fuscus) and the Eastern red bat (LABO = Lasiurus
borealis), which are stored in the vignettes/plots/
folder:
????_state_fatalities_animation.gif
????_annual_fatalities.png
????_cumulative_fatalities.png
Finally, we can consider how relative risk for each species may vary across the U.S. by pairing our predictions with projections from the National Renewable Energy Laboratory (NREL) that describe how much wind energy could potentially be added to the landscape under different scenarios.
We predict the number of nationwide bat fatalities that could occur in association with varying scenarios of wind energy build-out using the code below:
#Determine the species name (four letter code) from the species_columns
species_names <- substr(species_columns, 1, 4)
#For each selected species, create relative risk maps of potential fatalities
for(species in seq(1,length(species_names))) {
#Select the species
species_name <- species_names[species]
#Print some messages to provide progress feedback to user
print(paste0("Processing species: ", species_name))
#Determine the range of the species
species_range <- spc_ranges |> filter(Spp_code == species_name)
#Load the center/scale information for the covariates
cs_save <- jags_data$cs_save
#Create a data.frame of parameter values to run predictions on
parameter_grid <- data.frame(intercept = 1,
prop_spring = 1,
prop_summer = 1,
prop_fall = 1,
cutin = 0,
max_cutin_speed = 3,
year = 2024)
#Create a subset of facilities added after 2020 representing new turbines
facilities_post_2020 <- facilities |>
filter(start_year >= 2020)
#Print a progress note
print("Intersecting NREL grid with species range... this could take 5+ minutes!")
#Determine which nrel_grid cells fall within species range
nrel_grid_in_range <- NationalBatsAndWind::nrel_grid |> st_filter(species_range)
#For each NREL scenario, create a plot of relative risk
for(scenario in c("limited", "reference", "open")) {
#Print another note
print(paste0("Creating map for ", scenario, " scenario"))
#Open the corresponding scenario grid
scenario_grid <- read.csv(paste0(
"../data-raw/grid/nrel_grid/",
scenario,
"_2030_moderate_115hh_170rd_supply-curve.csv")) |>
dplyr::mutate(fid = sc_point_gid)
#Determine covariate values for hypothetical facilities at each grid cell
nrel_turbines_summary <- nrel_grid_in_range |>
left_join(scenario_grid, by = "fid") |>
mutate(
year = parameter_grid$year[1],
age = 0, #what if we put a new facility at this location?
p_cap = capacity_mw,
sum_t_cap = capacity_mw,
mean_ttlh = hub_height + 170/2,
mean_hh = hub_height,
total_rsa = n_turbines * (170 / 2)^2 * 3.14,
mean_rsa = (170 / 2)^2 * 3.14,
n_turbine = n_turbines,
total_mw = p_cap,
years_since_wns = year - wns_yoa_av,
years_since_wns_nonneg = if_else(years_since_wns < 0, 0, years_since_wns),
mean_lrs = hub_height - 170/2,
lrs_mean_windspeed_av = windspeed_10_m_av + (2/3) * (windspeed_40_m_av - windspeed_10_m_av),
shear_ttlh_minus_lrs_av = windspeed_200_m_av - (windspeed_10_m_av + (2/3) * (windspeed_40_m_av - windspeed_10_m_av)),
total_rsa_pct = total_rsa,
shape_ind = median(facilities_post_2020$shape_ind),
para = median(facilities_post_2020$para),
mean_dist_to_edge = median(facilities_post_2020$mean_dist_to_edge),
intercept = parameter_grid[1,1],
prop_spring = parameter_grid[1,2],
prop_summer = parameter_grid[1,3],
prop_fall = parameter_grid[1,4],
vapor_pr_av = vapor_pr_av,
vapor_pr_av_sq = NA,
cutin = parameter_grid[1,5],
max_cutin_speed = parameter_grid[1,6],
physdiv_av = physdiv_av)
#Center and scale the covariates according to cs_save
for(i in which(rownames(cs_save) %in% colnames(jags_data$X))){
if(class(nrel_turbines_summary[[rownames(cs_save)[i]]])=="units"){
nrel_turbines_summary[[rownames(cs_save)[i]]] <- units::drop_units(nrel_turbines_summary[[rownames(cs_save)[i]]])
}
nrel_turbines_summary <- nrel_turbines_summary %>%
mutate(!!rownames(cs_save)[i] := (!!sym(rownames(cs_save)[i]) - cs_save[i,1] )/cs_save[i,2])
}
#After center/scaling, calculate the quadratic effects for a few covariates
nrel_turbines_summary <- nrel_turbines_summary |>
mutate(vapor_pr_av_sq = vapor_pr_av^2)
#Predict across the new cells
temp <- nrel_turbines_summary %>%
NationalBatsAndWind::multispecies_predict(
jags_data,
jags_fit,
species = species,
fast = T)
#Calculate the 0.025, 0.5, and 0.975 quantiles and mean of the posteriors
#across all fid cells included in the species range from the nrel_grid
predno_q2.5 <- paste0('pred_q2.5_', parameter_grid$year[1])
predno_q50 <- paste0('pred_q50_', parameter_grid$year[1])
predno_q97.5 <- paste0('pred_q97.5_', parameter_grid$year[1])
out <- data.frame(fid = temp$cells_modeled)
out <- out %>%
mutate(!!predno_q2.5 := apply(temp$pred, 1, quantile, 0.025, na.rm = TRUE)) %>%
mutate(!!predno_q50 := apply(temp$pred, 1, quantile, 0.5, na.rm = TRUE)) %>%
mutate(!!predno_q97.5 := apply(temp$pred, 1, quantile, 0.975, na.rm = TRUE))
#Add predictions to a new copy of the nrel_grid
nrel_grid_pred <- nrel_turbines_summary
nrel_grid_pred <- nrel_grid_pred |>
left_join(out, by = 'fid')
#Create a map of the relative risk of potential fatalities for the species
potential_fatalities <- ggplot() +
geom_sf(
data = nrel_grid_pred,
aes(fill = if_else(
pred_q50_2024 > 5000,
5000,
if_else(pred_q50_2024 < 1, 1, pred_q50_2024))),
color = "transparent") +
scale_fill_gradientn(
trans = "log",
colors = gist_ncar(100),
labels = scales::unit_format(unit = "", big.mark = ",", accuracy = NULL),
na.value = "gray25",
breaks = c(5,50,500,5000),
limits = c(1,5000)) +
geom_polygon(
data = states,
aes(x = x, y = y, group = group),
color = "black",
fill = "transparent",
linewidth = 0.25) +
guides(shape = guide_legend(
"Fatality report",
override.aes = list(shape=16, size = 2))) +
guides(size = "none") +
ggthemes::theme_map() +
labs(fill = paste0("Annual ", species_name, "\nfatality potential")) +
labs(size = NULL) +
theme_minimal() +
theme(
plot.background = element_rect(fill = "white", color = "white"),
legend.position = "right",
text = element_text(size=16, colour="black"),
axis.text = element_text(size=12, colour="black"),
axis.title = element_text(size=12,face="bold")) +
ylab("") +
xlab("") +
coord_sf(crs = crs(NationalBatsAndWind::nrel_grid))
#Save the map for later
ggsave(
paste0(
"plots/",
species_name, "_",
scenario, "_potential_bat_fatalities.png"),
plot = potential_fatalities,
width = 9,
height = 4,
units = c("in"),
dpi = 300,
create.dir = T
)
}
}
Again, let’s consider only the results for hoary bats in this vignette.
Under the “limited” NREL scenario, wind energy build-out is restricted by a high number of social, ecological, and economic constraints. Here, we plot hoary bat fatality potential across the United States. Fatality potential is calculated by predicting the potential maximum number of hoary bat fatalities within a specific grid cell, here, under the “limited” NREL scenario. Grayed out cells indicate that wind energy is restricted under this scenario and cannot be built in this location.
Under the “reference” NREL scenario, wind energy build-out is restricted by a moderate number of social, ecological, and economic constraints. Here, we plot hoary bat fatality potential for each grid cell under the “reference” NREL scenario.
Under the “open” NREL scenario, wind energy build-out has the lowest number of social, ecological, and economic constraints. Finally, we plot hoary bat fatality potential under the maximum amount of wind energy that might be built within each grid cell according to the “open” NREL scenario.
These maps meant to provide a general idea of where hoary bats may be at higher risk of fatalities due to wind energy build-out across the U.S.
Please review the plots that we have created for the other species,
big brown bat (EPFU) and Eastern red bat (LABO). Check the following
figures in the vignettes/plots/ folder:
????_limited_potential_bat_fatalities.png
????_reference_potential_bat_fatalities.png
????_open_potential_bat_fatalities.png
Here is the map for the “open” NREL scenario and annual potential big brown bat (EPFU) fatalities:
Finally, here is the map for the “open” NREL scenario and annual potential Eastern red bat (LABO) fatalities:
Again, please recall that this model evaluates only the publicly available data included in this release. In the near future we plan to publish our results on the full data set, which includes business confidential reports that we cannot share as raw data. Our analyses have shown that the results on the full data set are more robust and don’t necessarily agree fully with the example results we describe here. Additionally, we explore more covariate relationships, which refine our spatial predictions.